#include "includesDPC.h"
#include "pocsDPC.h"
// Forward projection projector.
using namespace std;
void P4fproject(ScanData & s)
{
	int i,j, N=-1;
	double x=0,total=0;
	vector<double> norm (s.Ny, 0);
	vector<double> T (s.Ny, 0);
	vector<double> F (s.Ny, 0);
	vector<double> D (s.Ny, 0);

	for(j=0 ; j < s.Ny ; j++){
		if(s.Lmask[j]){
			N++;
			for (i=0;i<s.Nx;i++){
				if (!s.APmask[i][j] && s.CSmask[i][j]) {
					T[N] = T[N] + s.dx*s.recon[i][j];
					x = s.xmin + i*s.dx;
					norm[N] = norm[N] + s.dx;
				}
			}
		}
	}
	
	// Differentiate
	for(j=0; j<N;++j){ 
		F[j] = (T[j+1] - T[j])/s.dy;
	}
	F[N] = F[N-1];

	N=-1;
	for(j = 0 ; j<s.Ny ; ++j){
		if(s.Lmask[j]){
			N++;
			T[N] = F[N];
			F[N] = (s.FPdata[j] - F[N])/norm[N];
		}
	}

	// Perform modified Euler
	N=-1;
	D[0]=0;
	for(j = 0; j<s.Ny-1; ++j){
		if(s.Lmask[j]){
			N++;
			D[N+1] = D[N] + s.dy*F[N];
			total = total + D[N+1];
		} 
	}
	//cout << total<<endl;
	N=-1;
	for(j=0; j<s.Ny ; j++){
		if(s.Lmask[j]){
			N++;
			//cout << D[N] - total/(N+1.)  <<", " << s.FPdata[j] << ", " << T[N] << ", "<<endl;
			T[N]=0;
			for (i=0 ; i < s.Nx ; i++){
				x = s.xmin + i*s.dx;
				if (!s.APmask[i][j] && s.CSmask[i][j]){
					s.recon[i][j] = s.recon[i][j] + D[N] - total/(N+1.);
				}
			}
		}
	}
	
}
